Objective

This vignette describes how copy number variation (CNV) can
be identified using Mosaic. Data filtering, normalization and
ploidy computation are described.

The h5 file used in this notebook can be found here

In [1]:
import missionbio.mosaic.io as mio

sample = mio.load('./sample.analyzed.h5')

Loading, ./sample.analyzed.h5
Loaded in 0.1s.

This is an analyzed h5 file. Hence the clones and clusters
are already labeled. It contains three cell lines KG-1,
Tom-1, and Jurkat. It also contains doublets of each pair
as seen in the heatmap.

For your h5 file, identify the clones using the DNA variants
before proceeding with this vignette.

In [2]:
sample.dna.heatmap('NGT')

Filtering amplicons

We will only be working with the DNA and CNV assays.

The CNV object contains 138 amplicons. But not all
of them worked as expected. We will drop amplicons
which worked in less than half the total cells.

In [3]:
# This returns the reads in each cell and amplicon

reads = sample.cnv.get_attribute('read_counts', constraint='row+col')
reads
Out[3]:
AMLV2_AMP1 AMLV2_AMP2 AMLV2_AMP3 AMLV2_AMP4 AMLV2_AMP5 AMLV2_AMP6 AMLV2_AMP7 AMLV2_AMP8 AMLV2_AMP9 AMLV2_AMP10 ... chr10_106721610 chr10_5554293 chr10_77210191 chr14_56969005 chr16_55770629 chr16_8569820 chr18_9750662 chr6_17076840 chr6_40116264 chr6_62094287
read_counts
AACAACCTAAACTTGTCG 139 161 105 219 92 157 152 127 74 178 ... 478 169 184 226 201 123 332 127 304 218
AACAACTGGTACGTTGGA 82 64 66 87 9 105 74 53 43 94 ... 255 132 100 134 124 31 155 126 173 102
AACAATGCAAGACCACGC 110 130 54 131 24 47 70 28 52 99 ... 419 32 37 181 139 98 170 322 306 102
AACAATGCACTGACGGAT 124 125 93 163 25 148 164 93 65 177 ... 495 58 22 124 233 78 195 197 307 214
AACAATGCAGAGCCATCC 50 59 31 39 14 51 49 16 29 82 ... 122 18 71 63 51 42 151 134 114 118
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
TTGTATCACGCTGTTCGG 38 17 16 5 3 67 22 11 12 27 ... 75 28 63 14 69 7 80 99 72 29
TTGTCAACCGCTTGTCCT 104 71 48 82 0 61 64 30 17 60 ... 187 25 127 86 96 48 74 134 153 104
TTGTCAACCTACAACACC 162 88 22 39 14 128 44 22 36 71 ... 379 53 132 122 235 125 294 153 345 82
TTGTCAACCTAGTAACGG 140 107 42 157 52 149 89 46 34 115 ... 189 96 142 177 153 120 265 207 274 89
TTGTTAGAGATCAGGATG 60 57 39 76 13 32 69 27 35 57 ... 267 10 87 55 104 82 95 111 151 96

2476 rows × 138 columns

In [4]:
# Only amplicons found in more than half the cells are analyzed
# The other amplicons are dropped.

working_amplicons = (reads.median() > 0).values
sample.cnv = sample.cnv[:, working_amplicons]
In [5]:
# We are now left with 135 amplicons

sample.cnv.shape
Out[5]:
(2476, 135)

Normalization and ploidy computation

Finding the ploidy of the cells can be split into
two distinct calculations.

  • Normalization

    This step is required to correct for systemic artefacts.
    Cell to cell, and amplicon to amplicon variations are corrected.
    Each cell is normalized to it's total reads and each amplicon
    to its median counts. This is performed using the normalize_reads
    method, which adds the normalized_counts layer to the CNV object.

  • Correcting the baseline

    Not all amplicons show copy loss or copy gain. In the normalization
    step all amplicons in a cell were normalized using the same value
    i.e. the total reads in the cell. This results in slight deviations
    from the true ploidy value of the cell for each amplicon.

    To correct for this, a clone which is known to be diploid is identified
    and all the other cells are normalized to the counts of this clone. This
    step is done using the compute_ploidy method, and providing it
    with the barcodes known to be diploid. This adds the ploidy layer
    to the CNV object.

In [6]:
# Reads are normalized to correct for systemic artefacts

sample.cnv.normalize_reads()
In [7]:
# Assuming Jurkat cells are diploid for all the amplicons,
# we can compute the ploidy of the other cell lines as follows

sample.cnv.compute_ploidy(diploid_cells=sample.dna.barcodes('Jurkat'))

Visualizations

The ploidy values are stored in the ploidy layer

These can be visualized using a heatmap or a ploidy line plot for each clone.

In [8]:
# Assign the DNA labels to the CNV object
# We want to ensure the labels are the same

sample.cnv.set_labels(sample.dna.get_labels())
sample.cnv.set_palette(sample.dna.get_palette())
In [9]:
# Heatmap with the features ordered by the default amplicon order

sample.cnv.heatmap('ploidy', features=sample.cnv.ids())

As seen in the heatmap, KG-1 has a copy loss for some
of the amplicons in the panel. TOM-1 does not show much
deviation from ploidy 2.

This can also be visualized by the per amplicon ploidy plot.
This plot is generated for each clone separately.

In [10]:
# KG-1 shows copy loss at two different locations

sample.cnv.plot_ploidy('KG-1')
In [11]:
# TOM-1 is diploid for all the amplicons

sample.cnv.plot_ploidy('TOM-1')
In [12]:
# Since Jurkat was used to set the baseline,
# it is perfectly diploid for all amplicons

sample.cnv.plot_ploidy('Jurkat')

Further analysis

Often copy loss or copy gain will be seen in a set of amplicons.
These amplicons might belong to the same gene. Cross checking the
amplicon positions with the gene positions can be useful.

Multiple amplicons showing a CNV makes the occurence much more
likely to be a true positive observation.